# Description ------------------------------------------------------------------
#   Analysis of zip codes common to both utilities
#
#   We want to consider three dimensions:
#     1. Type of price: 5 types
#     2. Number of lags: 0–3
#     3. Spatiotemporal fixed effects

# Setup ------------------------------------------------------------------------
  # Options
  options(stringsAsFactors = FALSE)
  # Packages
  library(pacman)
  p_load(stringr, magrittr)
  # Set the number of threads for the lfe package
  # options(lfe.threads = 32)
  # Directory of joined bill-weather files
  # dir_project <- "S:/Edward/Elasticity/"
  dir_project  <- "/Users/edwardarubin/Dropbox/Research/MyProjects/NaturalGas/"
  dir_rds      <- paste0(dir_project, "DataR/")
  dir_output   <- paste0(dir_project, "Text/ForSubmission/Tables/")
  dir_results  <- paste0(dir_rds, "Results/Results20171024/")
  dir_results2 <- paste0(dir_rds, "Results/Results20171025/")

# Functions for formatting rows ------------------------------------------------
  # Quick and simple formatting functions
  cf <- . %>% round(4) %>% format(nsmall = 4) %>% str_replace_all("-", "$-$")
  se <- . %>% format(trim = F, digits = 3) %>% paste0("(", ., ")")
  n <- . %>% format(trim = F, digits = 2, nsmall = 1,
    scientific = F, big.mark = ",")
  tf <- . %>% format(justify = "r")
  # Function to add LaTeX separators and endings
  tex_sep <- function(object) {
    # Number of columns in
    n_cols <- ncol(object)
    # Add separators
    lapply(
      X = 1:n_cols,
      FUN = function(j) {
        cbind(matrix(object[,j]), ifelse(j != n_cols, " & ", " \\\\ "))
      }) %>% do.call(cbind, .)
  }
  # Function to add significance stars (10%, 5%, 1%)
  add_stars <- function(coefs, p_vals) {
    # Add a star for each significance level
    the_stars <- paste0(
      ifelse(p_vals < 0.10, "*", ""),
      ifelse(p_vals < 0.05, "*", ""),
      ifelse(p_vals < 0.01, "*", ""))
    # Paste stars onto coefficients
    paste0(coefs, "$^{", the_stars, "}$")
  }

# Function: Simple, general tables ---------------------------------------------
  simple_table <- function(path = dir_output, filename, label,
    the_mat, row_names, row_ends = NULL, sideways = F, landscape = F,
    title_bf, title_nb, top_more = NULL, above_label,
    col_numbers = T, col_labels = NULL, the_caption = NULL) {
      # The dimensions of the matrix
      n_f <- ncol(the_mat)
      # Create table top
      the_top <- c(
        paste0("\\begin{", ifelse(sideways, "sideways", ""),
          "table}[!htbp] \\centering"),
        paste0("  \\caption{\\textbf{", title_bf, "}: ", title_nb, "}"),
        paste0("  \\label{", label, "}"),
        paste(c("\\begin{tabular}{@{\\extracolsep{3pt}}l",
          rep("c", n_f), "@{}}"), collapse = ""),
        switch(as.character(!is.null(above_label)),
          "TRUE" = paste0("  \\multicolumn{", n_f + 1, "}{c}{", above_label, "} \\\\"),
          "FALSE" = NULL),
        "  \\toprule\\midrule",
        "  \\ra{1.25}",
        NULL
      )
      # Add more (before column numbers/labels), if desired
      if (!is.null(top_more)) the_top %<>% c(top_more)
      # Column numbers, if desired
      if (col_numbers == T) {
        the_top <- paste0("  & ",
          paste0("(", 1:n_f, ")", collapse = " & "),
          " \\\\") %>% c(the_top, .)
      }
      # Add column labels, if desired
      if (!is.null(col_labels)) {
        the_top %<>% c(.,
          paste0(paste(col_labels, collapse = " & "), " \\\\ "))
      }
      # Mid rule
      the_top %<>% c(., "  \\midrule \\\\[-0.9em]")
      # Create the meat of the table
      the_meat <- cbind(matrix(row_names), the_mat) %>%
        # Spacing
        tex_sep() %>%
        # Add row ends
        cbind(., matrix(row_ends)) %>%
        # Formatting
        knitr::kable(format = "pandoc") %>%
        as.character() %>%
        tail(-1) %>% head(-1)
      # Create the bottom part of the table code
      the_bottom <- c(
        "\\bottomrule",
        "\\\\[-0.9em]",
        "\\end{tabular}",
        switch(as.character(is.null(the_caption)),
          "TRUE" = the_caption,
          "FALSE" = paste0(
            "\\caption*{\\raggedright\\footnotesize \\textbf{Notes:} ",
            the_caption, "}")),
        paste0("\\end{", ifelse(sideways, "sideways", ""), "table}"),
        NULL)
      # Put everything together
      the_out <- paste0(c(the_top, the_meat, the_bottom), "\n")
      # Add landscape command if "landscape" is requested
      if (landscape) {
        the_out <- c("\\begin{landscape}\n", the_out, "\\end{landscape}")
      }
      the_out <- capture.output(expr = cat(the_out))
      # Write it
      write(the_out, paste0(path, filename))
  }

# Function: First-stage tables -------------------------------------------------
  table_stage1 <- function(the_files, filename, path, sideways = F, landscape = F,
    title_bf, title_nb, label, top_label, col_labels, col_numbers = T,
    the_caption, hdd = T) {
    # Length of 'the_files' (columns)
    n_f <- length(the_files)
    # Create the top part of the table code
    the_top <- c(
      paste0("\\begin{", ifelse(sideways, "sideways", ""),
        "table}[!htbp] \\centering"),
      paste0("  \\caption{\\textbf{", title_bf, "}: ", title_nb, "}"),
      paste0("  \\label{", label, "}"),
      paste(c("\\begin{tabular}{@{\\extracolsep{3pt}}l",
        rep("c", n_f), "@{}}"), collapse = ""),
      "  \\toprule\\midrule",
      "  \\ra{1.25}",
      NULL
    )
    # Top label
    if (!is.null(top_label)) {
      the_top <- c(the_top,
        paste0("\\\\[-0.9em] & \\multicolumn{", n_f,
        "}{c}{\\textbf{", top_label, "}} \\\\ \\cmidrule{2-",
        n_f + 1, "}"))
    }
    # Column numbers, if desired
    if (col_numbers == T) {
      the_top <- paste0("  & ",
        paste0("(", seq_along(the_files), ")", collapse = " & "),
        " \\\\") %>% c(the_top, .)
    }
    # Add column labels, if desired
    if (!is.null(col_labels)) {
      the_top %<>% c(., paste0(paste(col_labels, collapse = " & "), " \\\\ "))
    }
    # Mid rule
    the_top %<>% c(., "  \\midrule \\\\[-0.9em]")
    # Create the bottom part of the table code
    the_bottom <- c(
      "\\bottomrule",
      "\\\\[-0.9em]",
      "\\end{tabular}",
      paste0("\\caption*{\\raggedright\\footnotesize \\textbf{Notes:} ",
        the_caption, "}"),
      paste0("\\end{", ifelse(sideways, "sideways", ""), "table}"),
      NULL
    )
    # Read results and create the body of the table
    tmp_results <- lapply(
      X = the_files,
      FUN = function(x) readRDS(paste0(dir_results, "/", x)))
    # Format and star first row of coefficients
    coef1 <- sapply(tmp_results, function(x) x %>% extract2(1) %$%coef %>%
      extract(2,1)) %>% cf()
    pval1 <- sapply(tmp_results, function(x) x %>% extract2(1) %$%coef %>%
      extract(2,4))
    coef1 %<>% add_stars(pval1)
    # Format and star second row of coefficients
    coef2 <- sapply(tmp_results, function(x) x %>% extract2(1) %$%coef %>%
      extract(3,1)) %>% cf()
    pval2 <- sapply(tmp_results, function(x) x %>% extract2(1) %$%coef %>%
      extract(3,4))
    coef2 %<>% add_stars(pval2)
    # Create the meat of the table
    the_meat <- rbind(
      # Coefficients from spot price
      coef1,
      # Standard errors from spot price
      sapply(tmp_results, function(x) x %>% extract2(1) %$%coef %>%
        extract(2,2)) %>% se(),
      # Coefficients from spot price interacted with utility
      coef2,
      # Standard errors from spot price interacted with utility
      sapply(tmp_results, function(x) x %>% extract2(1) %$%coef %>%
        extract(3,2)) %>% se(),
      # First-stage F statistic
      sapply(tmp_results, function(x) x %>% extract2(1) %$%f) %>% n(),
      # Bill HDD
      if (identical(hdd, T)) rep("T", n_f) %>% tf(),
      if (!identical(hdd, T)) hdd,
      # HH FE
      rep("T", n_f) %>% tf(),
      # City-month FE
      rep("T", n_f) %>% tf(),
      # Grab sample sizes
      sapply(tmp_results, function(x) x %>% extract2(1) %$%N) %>% n(),
      # NULL
      NULL
    ) %>%
    # Row names
    cbind(c("Spot price", "", "Spot price $\\times$ SoCalGas", "",
      "First-stage F stat.", "Bill HDDs", "Household FE", "City month-of-sample FE",
      "\\textit{N}"), .) %>%
    # LaTeX separators
    tex_sep() %>%
    # Print to screen nicely
    knitr::kable(format = "pandoc") %>% as.character() %>% tail(-1) %>% head(-1)
    # Add extra lines and spacing
    the_meat[2] <- paste0(the_meat[2], "\\\\[-0.5em]")
    the_meat[4] <- paste0(the_meat[4], "\\\\[-0.7em] \\midrule")
    the_meat %<>% paste0("  ", .)
    # Put everything together
    the_out <- paste0(c(the_top, the_meat, the_bottom), "\n")
    # Add landscape command if "landscape" is requested
    if (landscape) {
      the_out <- c("\\begin{landscape}\n", the_out, "\\end{landscape}")
    }
    the_out <- capture.output(expr = cat(the_out))
    # Write it
    write(the_out, paste0(path, filename))
  }

# Function: Second-stage tables ------------------------------------------------
  table_stage2 <- function(the_files, filename, path, sideways = F, landscape = F,
    title_bf, title_nb, label, above_label, top_label, top_custom = NULL,
    col_labels, col_labels2 = NULL, col_numbers = T, hdd = T, the_caption) {
    # Length of 'the_files' (columns)
    n_f <- length(the_files)
    # Create the top part of the table code
    the_top <- c(
      paste0("\\begin{", ifelse(sideways, "sideways", ""),
        "table}[!htbp] \\centering"),
      paste0("  \\caption{\\textbf{", title_bf, "}: ", title_nb, "}"),
      paste0("  \\label{", label, "}"),
      paste(c("\\begin{tabular}{@{\\extracolsep{3pt}}l",
        rep("c", n_f), "@{}}"), collapse = ""),
      ifelse(!is.null(above_label),
        paste0("  \\multicolumn{", n_f + 1, "}{c}{", above_label, "} \\\\"),
        NULL),
      "  \\toprule\\midrule",
      "  \\ra{1.25}",
      NULL
    )
    # Top label
    if (!is.null(top_label)) {
      the_top <- c(the_top,
        paste0("\\\\[-0.9em] & \\multicolumn{", n_f,
        "}{c}{\\textbf{", top_label, "}} \\\\ \\cmidrule{2-",
        n_f + 1, "}"))
    }
    # 'Custom' (user-defined) top label
    if (!is.null(top_custom)) the_top <- c(the_top, top_custom)
    # Column numbers, if desired
    if (col_numbers == T) {
      the_top <- paste0("  & ",
        paste0("(", seq_along(the_files), ")", collapse = " & "),
        " \\\\") %>% c(the_top, .)
    }
    # Add column labels, if desired
    if (!is.null(col_labels)) {
      the_top %<>% c(., paste0(paste(col_labels, collapse = " & "), " \\\\ "))
    }
    # Add second set of column labels, if desired
    if (!is.null(col_labels2)) {
      the_top %<>% c(., paste0(paste(col_labels2, collapse = " & "), " \\\\ "))
    }
    # Mid rule
    the_top %<>% c(., "  \\midrule \\\\[-0.9em]")
    # Create the bottom part of the table code
    the_bottom <- c(
      "\\bottomrule",
      "\\\\[-0.9em]",
      "\\end{tabular}",
      paste0("\\caption*{\\raggedright\\footnotesize \\textbf{Notes:} ",
        the_caption, "}"),
      paste0("\\end{", ifelse(sideways, "sideways", ""), "table}"),
      NULL
    )
    # Read results and create the body of the table
    tmp_results <- lapply(
      X = the_files,
      FUN = function(x) readRDS(paste0(dir_results, "/", x)))
    # Format and star first row of coefficients
    coef1 <- sapply(tmp_results, function(x) x %>% extract2(2) %$%coef %>%
      extract(2,1)) %>% cf()
    pval1 <- sapply(tmp_results, function(x) x %>% extract2(2) %$%coef %>%
      extract(2,4))
    coef1 %<>% add_stars(pval1)
    # Create the meat of the table
    the_meat <- rbind(
      # Coefficients from spot price
      coef1,
      # Standard errors from spot price
      sapply(tmp_results, function(x) x %>% extract2(2) %$%coef %>%
        extract(2,2)) %>% se(),
      # First-stage F statistic
      sapply(tmp_results, function(x) x %>% extract2(1) %$%f) %>% n(),
      # Bill HDD
      if (identical(hdd, T)) rep("T", n_f) %>% tf(),
      if (!identical(hdd, T)) hdd,
      # HH FE
      rep("T", n_f) %>% tf(),
      # City-month FE
      rep("T", n_f) %>% tf(),
      # Grab sample sizes
      sapply(tmp_results, function(x) x %>% extract2(1) %$%N) %>% n(),
      # NULL
      NULL
    ) %>%
    # Row names
    cbind(c(
      "Log(Price)",
      "\\quad \\textit{instrumented}",
      "First-stage F stat.", "Bill HDDs", "Household FE", "City month-of-sample FE",
      "\\textit{N}"), .) %>%
    # LaTeX separators
    tex_sep() %>%
    # Print to screen nicely
    knitr::kable(format = "pandoc") %>% as.character() %>% tail(-1) %>% head(-1)
    # Add extra lines and spacing
    the_meat[2] <- paste0(the_meat[2], "\\\\[-0.7em] \\midrule")
    the_meat %<>% paste0("  ", .)
    # Put everything together
    the_out <- paste0(c(the_top, the_meat, the_bottom), "\n")
    # Add landscape command if "landscape" is requested
    if (landscape) {
      the_out <- c("\\begin{landscape}\n", the_out, "\\end{landscape}")
    }
    the_out <- capture.output(expr = cat(the_out))
    # Write it
    write(the_out, paste0(path, filename))
  }

# Captions for results tables --------------------------------------------------
  cap_lag2 <- "Each price in the table is the second lag of price, \\textit{i.e.}, the prices from two bills prior to the current bill."
  cap_hub <- "\\textit{(HH) Spot price} refers to the weekly average spot price for natural gas at Louisiana's Henry Hub in the week preceding the utility's price change."
  cap_std <- "Each column denotes a separate regression. Errors are two-way clustered within (1) household and (2) utility by climate-zone by billing-cycle (the level at which price varies). All regressions include heating degree days (HDDs) within the households' bill."
  cap_avg <- "\\textit{Avg.} or \\textit{average} price is the total bill divided by quantity."
  cap_avgmrg <- "\\textit{Avg. Mrg.} or \\textit{average marginal} price denotes the quantity-weighted average of the household's marginal price."
  cap_sim <- "\\textit{Sim. Mrg.} or \\textit{simulated marginal} price is the household's marginal price (using the relevant pricing regime) as a function of the household's historical consumption patterns (lagged bills 10 through 14)."
  cap_base <- "\\textit{Base} or \\textit{baseline} price refers to the price the household pays for its first unit (\\textit{therm}) of natural gas."
  cap_nobs <- "As discussed in the empirical strategy section, the numbers of observations differ due to the lags required to calculate the \\textit{simulated instrument} for marginal price."
  cap_sig <- "\\textit{Significance levels:} *10\\%, **5\\%, ***1\\%."
  cap_lag <- "With regard to lags: \\textit{No lag} refers to the price for the household's contemporaneous bill; \\textit{1 Lag} refers to the price on the household's previous bill; \\textit{etc.}"
  cap_season <- "\\textit{Summer} includes April through September. \\textit{Winter} includes October through March."
  cap_care <- "\\textit{CARE} households participate in the California Alternative Rates for Energy (CARE) program. CARE targets low-income households and provides a 20 percent discount on natural gas bills."
  cap_het <- "We estimate the heterogeneity results by splitting the sample along the dimension(s) of heterogeneity and then individaully estimating the models."
  cap_nbr <- "\\textit{Common zips} refers the set of zip codes in which each zip code receives natural gas from both PG\\&E and SoCalGas. \\textit{Neighbors 1} includes the \\textit{common zips} and the zip codes that immediately neighbor the common zips. \\textit{Neighbors 2} adds the neighbors of these neighbors (adding the neighbors of \\textit{Neighbors 1}). \\textit{Neighbors 3} adds the neighbors of \\textit{Neighbors 2}. Figure~\\ref{fig:neighbors} depicts these sets of zip codes."

# tab:stage1_spot: First stage, all prices, lag 2 ------------------------------
  table_stage1(
    the_files = c(
      "LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvgMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceBase__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogMrgSimIv1014__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    filename = "stage1_spot.tex",
    path = dir_output,
    title_bf = "First-stage results",
    title_nb = "Instrumenting consumers' prices with Henry Hub spot price",
    label = "tab:stage1_spot",
    sideways = F,
    top_label = "Dependent variable: Log price",
    col_labels = c("  \\textbf{Type of Price:}",
      "Marginal", "Average", "Avg. Mrg.", "Baseline", "Sim. Mrg."),
    the_caption = paste(cap_hub, cap_avg, cap_avgmrg, cap_sim, cap_nobs, cap_std, cap_sig))

# tab:stage2_spot: Second stage, all prices, lag 2 -----------------------------
  table_stage2(
    the_files = c(
      "LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvgMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceBase__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogMrgSimIv1014__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot.tex",
    path = dir_output,
    title_bf = "Second-stage results",
    title_nb = "Instrumenting consumers' prices with Henry Hub spot price",
    label = "tab:stage2_spot",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Type of Price",
    col_labels = c("  ",
      "Marginal", "Average", "Avg. Mrg.", "Baseline", "Sim. Mrg."),
    the_caption = paste(cap_lag2, cap_avg, cap_avgmrg, cap_sim, cap_nobs, cap_hub, cap_std, cap_sig))

# tab:stage12_spot: Both stages, all prices, lag 2 -----------------------------
  # Load data (both stages)
  tmp_mat <- lapply(X = c(
      "LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvgMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceBase__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogMrgSimIv1014__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # First stage: Grab and format the coeficient and standard error
      the_cf1 <- tmp[[1]]$coef[2:3,1] %>% cf()
      the_se1 <- tmp[[1]]$coef[2:3,2] %>% se()
      the_cf1 %<>% add_stars(tmp[[1]]$coef[2:3,4])
      # Second stage: Grab and format the coeficient and standard error
      the_cf2 <- tmp[[2]]$coef[2,1] %>% cf()
      the_se2 <- tmp[[2]]$coef[2,2] %>% se()
      the_cf2 %<>% add_stars(tmp[[2]]$coef[2,4])
      # Sample size and F stat.
      the_f <- tmp[[1]]$f %>% n()
      the_n <- tmp[[2]]$N %>% n()
      # Return a column matrix of the results
      c(the_cf1[1], the_se1[1], the_cf1[2], the_se1[2],
        the_cf2, the_se2,
        the_f, rep("T", 3), the_n) %>% matrix(ncol = 1)
    }) %>% do.call("cbind", .)
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage12_spot.tex",
    label = "tab:stage12_spot",
    the_mat = tmp_mat,
    row_names = c(
      "\\quad Spot price",
      "",
      "\\quad Spot price",
      "\\quad\\quad$\\times$ SoCalGas",
      "\\quad \\text{Log(Price)}",
      "\\quad\\quad \\textit{(instrumented)}",
      "First-stage F stat.",
      "Bill HDDs",
      "Household FE",
      "City mo.-of-sample FE",
      "\\textit{N}",
      NULL),
    row_ends = c(
      "",
      "\\\\[-0.5em]",
      "",
      "\\\\[-0.9em] \\midrule \n \\multicolumn{6}{l}{\\hspace{-1ex}\\textbf{Panel B:} Second-stage results} \\\\[0.5em]",
      "",
      "\\\\[-0.9em] \\midrule",
      rep("", 5)),
    sideways = F,
    title_bf = "First- and second-stage results",
    title_nb = "\\newline Instrumenting consumers' prices with the Henry Hub spot price",
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_more = "\\\\[-1em] \n \\multicolumn{6}{l}{\\hspace{-1ex}\\textbf{Panel A:} First-stage results} \\\\",
    col_numbers = T,
    col_labels = c("", "Marginal", "Average", "Avg. Mrg.", "Baseline", "Sim. Mrg."),
    the_caption = paste(cap_std, cap_hub, cap_lag2, cap_avg, cap_avgmrg, cap_sim, cap_nobs, cap_sig))

# tab:stage2_spot_mrg: Second stage, mrg price, lags -1-3 ----------------------
  table_stage2(
    the_files = c(
      "LogPriceMrg__Lead1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceMrg__Lag0__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceMrg__Lag1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceMrg__Lag3__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_mrg.tex",
    path = dir_output,
    title_bf = "Comparing lags, second-stage results",
    title_nb = "Marginal price with HH spot price IV",
    label = "tab:stage2_spot_mrg",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Lag of Marginal Price",
    col_labels = c("  ",
      "1 Lead", "No lag", "1 Lag", "2 Lags", "3 Lags"),
    the_caption = paste(cap_lag, cap_hub, cap_std, cap_sig))

# tab:stage2_spot_sim: Second stage, sim mrg price, lags -1-3 ------------------
  table_stage2(
    the_files = c(
      "LogMrgSimIv1014__Lead1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogMrgSimIv1014__Lag0__CityYm_HhId__All__All__BillHddDays.rds",
      "LogMrgSimIv1014__Lag1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogMrgSimIv1014__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogMrgSimIv1014__Lag3__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_sim.tex",
    path = dir_output,
    title_bf = "Comparing lags, second-stage results",
    title_nb = "Sim. marginal price with HH spot price IV",
    label = "tab:stage2_spot_sim",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Lag of Simulated Marginal Price",
    col_labels = c("  ",
      "1 Lead", "No lag", "1 Lag", "2 Lags", "3 Lags"),
    the_caption = paste(cap_lag, cap_sim, cap_hub, cap_std, cap_sig))

# tab:stage2_spot_avgmrg: Second stage, avg mrg price, lags -1-3 ---------------
  table_stage2(
    the_files = c(
      "LogPriceAvgMrg__Lead1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvgMrg__Lag0__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvgMrg__Lag1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvgMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvgMrg__Lag3__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_avgmrg.tex",
    path = dir_output,
    title_bf = "Comparing lags, second-stage results",
    title_nb = "Avg. marginal price with HH spot price IV",
    label = "tab:stage2_spot_avgmrg",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Lag of Average Marginal Price",
    col_labels = c("  ",
      "1 Lead", "No lag", "1 Lag", "2 Lags", "3 Lags"),
    the_caption = paste(cap_lag, cap_sim, cap_hub, cap_std, cap_sig))

# tab:stage2_spot_avg: Second stage, avg price, lags -1-3 ----------------------
  table_stage2(
    the_files = c(
      "LogPriceAvg__Lead1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag0__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag3__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_avg.tex",
    path = dir_output,
    title_bf = "Comparing lags, second-stage results",
    title_nb = "Avgerage price with HH spot price IV",
    label = "tab:stage2_spot_avg",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Lag of Average Price",
    col_labels = c("  ",
      "1 Lead", "No lag", "1 Lag", "2 Lags", "3 Lags"),
    the_caption = paste(cap_lag, cap_avg, cap_hub, cap_std, cap_sig))

# tab:stage2_spot_base: Second stage, baseline price, lags -1-3 ----------------
  table_stage2(
    the_files = c(
      "LogPriceBase__Lead1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceBase__Lag0__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceBase__Lag1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceBase__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceBase__Lag3__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_base.tex",
    path = dir_output,
    title_bf = "Comparing lags, second-stage results",
    title_nb = "Baseline price with HH spot price IV",
    label = "tab:stage2_spot_base",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Lag of Baseline Price",
    col_labels = c("  ",
      "1 Lead", "No lag", "1 Lag", "2 Lags", "3 Lags"),
    the_caption = paste(cap_lag, cap_base, cap_hub, cap_std, cap_sig))

# tab:stage2_spot_mrg_avg: Second stage, mrg and avg, lags -1-2 -----------------
  table_stage2(
    the_files = c(
      "LogPriceMrg__Lead1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceMrg__Lag0__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceMrg__Lag1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lead1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag0__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag1__CityYm_HhId__All__All__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_mrg_avg.tex",
    path = dir_output,
    title_bf = "Comparing lags, second-stage results",
    title_nb = "Marginal and average prices with HH spot price IV",
    label = "tab:stage2_spot_mrg_avg",
    sideways = F,
    landscape = T,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = NULL,
    top_custom = paste("  &",
      "\\multicolumn{4}{c}{\\textbf{Marginal Price}} &",
      "\\multicolumn{4}{c}{\\textbf{Average Price}} \\\\",
      "\\cmidrule(lr{1pt}){2-5} \\cmidrule(lr{1pt}){6-9}"),
    col_labels = c("  ",
      "1 Lead", "No lag", "1 Lag", "2 Lags", "1 Lead", "No lag", "1 Lag", "2 Lags"),
    the_caption = paste(cap_lag, cap_avg, cap_hub, cap_std, cap_sig))

# tab:stage2_spot_mrg_het1: 1-dimensional het, mrg, lag 2 -----------------------
  table_stage2(
    the_files = c(
      "LogPriceMrg__Lag2__CityYm_HhId__All__Summer__BillHddDays.rds",
      "LogPriceMrg__Lag2__CityYm_HhId__All__Winter__BillHddDays.rds",
      "LogPriceMrg__Lag2__CityYm_HhId__Care__All__BillHddDays.rds",
      "LogPriceMrg__Lag2__CityYm_HhId__Noncare__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_mrg_het1.tex",
    path = dir_output,
    title_bf = "Heterogeneity by season or income",
    title_nb = "\\newline Second-stage results, instrumenting marginal price with HH spot price",
    label = "tab:stage2_spot_mrg_het1",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Marginal Price",
    top_custom = paste("  &",
      "\\multicolumn{2}{c}{Split by \\textbf{Season}} &",
      "\\multicolumn{2}{c}{Split by \\textbf{CARE (Income)}} \\\\",
      "\\cmidrule(lr{1pt}){2-3} \\cmidrule(lr{1pt}){4-5}"),
    col_labels = c("  ",
      "Summer", "Winter", "CARE", "Non-CARE"),
    the_caption = paste(cap_std, cap_lag2, cap_season, cap_care, cap_het, cap_sig))

# tab:stage2_spot_avg_het1: 1-dimensional het, avg, lag 2 -----------------------
  table_stage2(
    the_files = c(
      "LogPriceAvg__Lag2__CityYm_HhId__All__Summer__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__All__Winter__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__Care__All__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__Noncare__All__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_avg_het1.tex",
    path = dir_output,
    title_bf = "Heterogeneity by season or income",
    title_nb = "\\newline Second-stage results, instrumenting average price with HH spot price",
    label = "tab:stage2_spot_avg_het1",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Average Price",
    top_custom = paste("  &",
      "\\multicolumn{2}{c}{Split by \\textbf{Season}} &",
      "\\multicolumn{2}{c}{Split by \\textbf{CARE (Income)}} \\\\",
      "\\cmidrule(lr{1pt}){2-3} \\cmidrule(lr{1pt}){4-5}"),
    col_labels = c("  ",
      "Summer", "Winter", "CARE", "Non-CARE"),
    the_caption = paste(cap_std, cap_season, cap_care, cap_het, cap_avg, cap_lag2, cap_sig))

# tab:stage2_spot_mrg_het2: 2-dimensional het, mrg, lag 2 -----------------------
  table_stage2(
    the_files = c(
      "LogPriceMrg__Lag2__CityYm_HhId__Care__Summer__BillHddDays.rds",
      "LogPriceMrg__Lag2__CityYm_HhId__Noncare__Summer__BillHddDays.rds",
      "LogPriceMrg__Lag2__CityYm_HhId__Care__Winter__BillHddDays.rds",
      "LogPriceMrg__Lag2__CityYm_HhId__Noncare__Winter__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_mrg_het2.tex",
    path = dir_output,
    title_bf = "Heterogeneity by season and income",
    title_nb = "\\newline Second-stage results, instrumenting marginal price with HH spot price",
    label = "tab:stage2_spot_mrg_het2",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Marginal Price",
    col_labels = c("  ",
      "Summer", "Summer", "Winter", "Winter"),
    col_labels2 = c("  ",
      "CARE", "Non-CARE", "CARE", "Non-CARE"),
    the_caption = paste(cap_std, cap_lag2, cap_season, cap_care, cap_het, cap_sig))

# tab:stage2_spot_avg_het2: 2-dimensional het, avg, lag 2 -----------------------
  table_stage2(
    the_files = c(
      "LogPriceAvg__Lag2__CityYm_HhId__Care__Summer__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__Noncare__Summer__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__Care__Winter__BillHddDays.rds",
      "LogPriceAvg__Lag2__CityYm_HhId__Noncare__Winter__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_avg_het2.tex",
    path = dir_output,
    title_bf = "Heterogeneity by season and income",
    title_nb = "\\newline Second-stage results, instrumenting average price with HH spot price",
    label = "tab:stage2_spot_avg_het2",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Average Price",
    col_labels = c("  ",
      "Summer", "Summer", "Winter", "Winter"),
    col_labels2 = c("  ",
      "CARE", "Non-CARE", "CARE", "Non-CARE"),
    the_caption = paste(cap_std, cap_lag2, cap_season, cap_care, cap_het, cap_avg, cap_lag2, cap_sig))

# tab:stage2_spot_sim_het2: 2-dimensional het, sim mrg, lag 2 -------------------
  table_stage2(
    the_files = c(
      "LogMrgSimIv1014__Lag2__CityYm_HhId__Care__Summer__BillHddDays.rds",
      "LogMrgSimIv1014__Lag2__CityYm_HhId__Noncare__Summer__BillHddDays.rds",
      "LogMrgSimIv1014__Lag2__CityYm_HhId__Care__Winter__BillHddDays.rds",
      "LogMrgSimIv1014__Lag2__CityYm_HhId__Noncare__Winter__BillHddDays.rds",
      NULL),
    filename = "stage2_spot_sim_het2.tex",
    path = dir_output,
    title_bf = "Heterogeneity by season and income",
    title_nb = "\\newline Second-stage results, instrumenting sim. mrg. price with HH spot price",
    label = "tab:stage2_spot_sim_het2",
    sideways = F,
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_label = "Simulated Marginal Price",
    col_labels = c("  ",
      "Summer", "Summer", "Winter", "Winter"),
    col_labels2 = c("  ",
      "CARE", "Non-CARE", "CARE", "Non-CARE"),
    the_caption = paste(cap_std, cap_lag2, cap_season, cap_care, cap_het, cap_sim, cap_lag2, cap_sig))

# tab:stage1_sim_mrg: Mrg price regressed on the sim mrg price  ----------------
  # Load the results and create the matrix for the table
  tmp_mat <- lapply(X = c("SIM__10_14__Lag0.rds", "SIM__11_13__Lag0.rds"),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[2]]$coef[1,1] %>% cf()
      the_se <- tmp[[2]]$coef[1,2] %>% se()
      the_cf %<>% add_stars(tmp[[2]]$coef[1,4])
      the_n <- tmp[[2]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf, the_se, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Add rows for description
  tmp_mat <- rbind(tmp_mat[1:2,],
    # HDDs
    c("T","T"),
    # Household FE
    c("T","T"),
    # City month-of-sample FE
    c("T","T"),
    # Lags
    c("10--14", "11--13"),
  tmp_mat[3,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage1_sim_mrg.tex",
    label = "tab:stage1_sim_mrg",
    the_mat = tmp_mat,
    row_names = c(
      "Simulated marginal price",
      "",
      "Bill HDDs",
      "Household FE",
      "City month-of-sample FE",
      "Lags used for sim. inst.",
      "\\textit{N}"),
    row_ends = c("", "\\\\[-0.7em] \\midrule", rep("", 5)),
    sideways = F,
    title_bf = "Testing the simulated instrument",
    title_nb = "\\newline Regressing marginal price on \\textit{simulated} marginal price",
    above_label = "\\textbf{Dependent variable:} Marginal price",
    top_more = NULL,
    col_numbers = T,
    col_labels = NULL,
    the_caption = paste(cap_std, str_replace(cap_sim, "14", "14 or 11 through 13"), cap_nobs, cap_sig))

# tab:ols: OLS results ---------------------------------------------------------
  # Load the results and create the matrix for the table
  tmp_mat <- lapply(X = c(
    "JointCA/jointCA__OLS__LogThmDay__LogPriceMrg_BillHddDays__Lag2__Ym_HhId__All__All.rds",
    "OLS__LogThmDay__LogPriceMrg_BillHddDays__Lag2__CityYm_HhId__All__All.rds",
    "OLS__LogThmDay__LogPriceMrg_BillHddDays__Lag2__Ym_HhId__All__All.rds",
    "JointCA/jointCA__OLS__LogThmDay__LogPriceBase_BillHddDays__Lag2__Ym_HhId__All__All.rds",
    "OLS__LogThmDay__LogPriceBase_BillHddDays__Lag2__CityYm_HhId__All__All.rds",
    "OLS__LogThmDay__LogPriceBase_BillHddDays__Lag2__Ym_HhId__All__All.rds",
    NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[2]]$coef[1,1] %>% cf()
      the_se <- tmp[[2]]$coef[1,2] %>% se()
      the_cf %<>% add_stars(tmp[[2]]$coef[1,4])
      the_n <- tmp[[2]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf, the_se, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Change table's structure: two prices get different rows
  tmp_mat <- rbind(
    cbind(tmp_mat[1:2,1:3], matrix(data = "", nrow = 2, ncol = 3)),
    cbind(matrix(data = "", nrow = 2, ncol = 3), tmp_mat[1:2,4:6]),
    tmp_mat[3,],
    NULL)
  # Add rows for description
  tmp_mat <- rbind(tmp_mat[1:4,],
    # HDDs
    c(rep("T", 6)),
    # Household FE
    c(rep("T", 6)),
    # Month-of-sample FE
    c("T","T","F","T","T","F"),
    # City month-of-sample FE
    c("F","F","T","F","F","T"),
    # Dataset
    c("5\\% CA", "Border", "Border", "5\\% CA", "Border", "Border"),
  tmp_mat[5,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "ols.tex",
    label = "tab:ols",
    the_mat = tmp_mat,
    row_names = c(
      "Log(Marginal price)",
      "",
      "Log(Baseline price)",
      "",
      "Bill HDDs",
      "Household FE",
      "Month-of-sample FE",
      "City by month-of-sample FE",
      "Sample",
      "\\textit{N}"),
    row_ends = c("", "\\\\[-0.7em]", "", "\\\\[-0.7em] \\midrule", rep("", 6)),
    sideways = T,
    title_bf = "OLS Results",
    title_nb = "Estimating elasticities, varying the dataset, price, and fixed effects",
    above_label = NULL,
    top_more = " & \\multicolumn{6}{c}{\\textbf{Dependent variable:} Log(Consumption, daily avg.)} \\\\ \\cmidrule{2-7}",
    col_numbers = T,
    col_labels = NULL,
    the_caption = paste(cap_std, cap_base, cap_lag2, cap_sig))

# tab:stage1_mrg_rob: First stage, mrg price, lag 2, vary HDDs and FEs ---------
  tmp_mat <- lapply(X = c(
    "LogPriceMrg__Lag2__CityYm_HhId__All__All__1.rds",
    "LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
    "LogPriceMrg__Lag2__CityYw_HhId__All__All__BillHddDays.rds",
    "LogPriceMrg__Lag2__ZipYw_HhId__All__All__BillHddDays.rds",
    NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # Determine whether we have HDDs
      nc <- grepl("hdd", a_file, ignore.case = T) + 1
      nc <- nc:(nc+1)
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[1]]$coef[nc,1] %>% cf()
      the_se <- tmp[[1]]$coef[nc,2] %>% se()
      the_cf %<>% add_stars(tmp[[1]]$coef[nc,4])
      the_f <- tmp[[1]]$f %>% n()
      the_n <- tmp[[1]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf[1], the_se[1], the_cf[2], the_se[2], the_f, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Add more rows
  tmp_mat <- rbind(
    tmp_mat[1:5,],
    # HDDs
    c("F", rep("T", 3)),
    # Household FE
    c("T", "T", "T", "T"),
    # City month-of-sample FE
    c("T", "T", "F", "F"),
    # City week-of-sample FE
    c("F", "F", "T", "F"),
    # Zip week-of-sample FE
    c("F", "F", "F", "T"),
    # N
    tmp_mat[6,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage1_mrg_rob.tex",
    label = "tab:stage1_mrg_rob",
    the_mat = tmp_mat,
    row_names = c(
      "Spot price",
      "",
      "Spot price $\\times$ SoCalGas",
      "",
      "First-stage F stat.",
      "Bill HDDs",
      "Household FE",
      "City by month-of-sample FE",
      "City by week-of-sample FE",
      "Zip by week-of-sample FE",
      "\\textit{N}"),
    row_ends = c(
      "",
      "\\\\[-0.7em]",
      "",
      "\\\\[-0.9em] \\midrule",
      rep("", 7)),
    sideways = F,
    title_bf = "First-stage results",
    title_nb = "\\newline Robustness to specification: Marginal price instrumented with spot price",
    above_label = "\\textbf{Dependent variable:} Log(Marginal price)",
    top_more = NULL,
    col_numbers = T,
    col_labels = NULL,
    the_caption = paste(cap_std, cap_lag2, cap_hub, cap_sig))

# tab:stage2_mrg_rob: Second stage, mrg price, lag 2, vary HDDs and FEs --------
  tmp_mat <- lapply(X = c(
    "LogPriceMrg__Lag2__CityYm_HhId__All__All__1.rds",
    "LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
    "LogPriceMrg__Lag2__CityYw_HhId__All__All__BillHddDays.rds",
    "LogPriceMrg__Lag2__ZipYw_HhId__All__All__BillHddDays.rds",
    NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # Determine whether we have HDDs
      nc <- grepl("hdd", a_file, ignore.case = T) + 1
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[2]]$coef[nc,1] %>% cf()
      the_se <- tmp[[2]]$coef[nc,2] %>% se()
      the_cf %<>% add_stars(tmp[[2]]$coef[nc,4])
      the_f <- tmp[[1]]$f %>% n()
      the_n <- tmp[[1]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf, the_se, the_f, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Add more rows
  tmp_mat <- rbind(
    tmp_mat[1:3,],
    # HDDs
    c("F", rep("T", 3)),
    # Household FE
    c("T", "T", "T", "T"),
    # City month-of-sample FE
    c("T", "T", "F", "F"),
    # City week-of-sample FE
    c("F", "F", "T", "F"),
    # Zip week-of-sample FE
    c("F", "F", "F", "T"),
    # N
    tmp_mat[4,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage2_mrg_rob.tex",
    label = "tab:stage2_mrg_rob",
    the_mat = tmp_mat,
    row_names = c(
      "Log(Marginal price)",
      "\\quad\\textit{instrumented}",
      "First-stage F stat.",
      "Bill HDDs",
      "Household FE",
      "City by month-of-sample FE",
      "City by week-of-sample FE",
      "Zip by week-of-sample FE",
      "\\textit{N}"),
    row_ends = c(
      "",
      "\\\\[-0.9em] \\midrule",
      rep("", 7)),
    sideways = F,
    title_bf = "Second-stage results",
    title_nb = "\\newline Robustness to specification: Marginal price instrumented with spot price",
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_more = NULL,
    col_numbers = T,
    col_labels = NULL,
    the_caption = paste(cap_std, cap_lag2, cap_sig))

# tab:stage2_avg_rob: Second stage, avg price, lag 2, vary HDDs and FEs --------
  tmp_mat <- lapply(X = c(
    "LogPriceAvg__Lag2__CityYm_HhId__All__All__1.rds",
    "LogPriceAvg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
    "LogPriceAvg__Lag2__CityYw_HhId__All__All__BillHddDays.rds",
    "LogPriceAvg__Lag2__ZipYw_HhId__All__All__BillHddDays.rds",
    NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # Determine whether we have HDDs
      nc <- grepl("hdd", a_file, ignore.case = T) + 1
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[2]]$coef[nc,1] %>% cf()
      the_se <- tmp[[2]]$coef[nc,2] %>% se()
      the_cf %<>% add_stars(tmp[[2]]$coef[nc,4])
      the_f <- tmp[[1]]$f %>% n()
      the_n <- tmp[[1]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf, the_se, the_f, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Add more rows
  tmp_mat <- rbind(
    tmp_mat[1:3,],
    # HDDs
    c("F", rep("T", 3)),
    # Household FE
    c("T", "T", "T", "T"),
    # City month-of-sample FE
    c("T", "T", "F", "F"),
    # City week-of-sample FE
    c("F", "F", "T", "F"),
    # Zip week-of-sample FE
    c("F", "F", "F", "T"),
    # N
    tmp_mat[4,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage2_avg_rob.tex",
    label = "tab:stage2_avg_rob",
    the_mat = tmp_mat,
    row_names = c(
      "Log(Average price)",
      "\\quad\\textit{instrumented}",
      "First-stage F stat.",
      "Bill HDDs",
      "Household FE",
      "City by month-of-sample FE",
      "City by week-of-sample FE",
      "Zip by week-of-sample FE",
      "\\textit{N}"),
    row_ends = c(
      "",
      "\\\\[-0.9em] \\midrule",
      rep("", 7)),
    sideways = F,
    title_bf = "Second-stage results",
    title_nb = "\\newline Robustness to specification: Average price instrumented with spot price",
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_more = NULL,
    col_numbers = T,
    col_labels = NULL,
    the_caption = paste(cap_std, cap_lag2, cap_avg, cap_sig))

# tab:stage2_avgmrg_rob: Second stage, avg price, lag 2, vary HDDs and FEs --------
  tmp_mat <- lapply(X = c(
    "LogPriceAvgMrg__Lag2__CityYm_HhId__All__All__1.rds",
    "LogPriceAvgMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
    "LogPriceAvgMrg__Lag2__CityYw_HhId__All__All__BillHddDays.rds",
    "LogPriceAvgMrg__Lag2__ZipYw_HhId__All__All__BillHddDays.rds",
    NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # Determine whether we have HDDs
      nc <- grepl("hdd", a_file, ignore.case = T) + 1
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[2]]$coef[nc,1] %>% cf()
      the_se <- tmp[[2]]$coef[nc,2] %>% se()
      the_cf %<>% add_stars(tmp[[2]]$coef[nc,4])
      the_f <- tmp[[1]]$f %>% n()
      the_n <- tmp[[1]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf, the_se, the_f, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Add more rows
  tmp_mat <- rbind(
    tmp_mat[1:3,],
    # HDDs
    c("F", rep("T", 3)),
    # Household FE
    c("T", "T", "T", "T"),
    # City month-of-sample FE
    c("T", "T", "F", "F"),
    # City week-of-sample FE
    c("F", "F", "T", "F"),
    # Zip week-of-sample FE
    c("F", "F", "F", "T"),
    # N
    tmp_mat[4,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage2_avgmrg_rob.tex",
    label = "tab:stage2_avgmrg_rob",
    the_mat = tmp_mat,
    row_names = c(
      "Log(Avg. marginal price)",
      "\\quad\\textit{instrumented}",
      "First-stage F stat.",
      "Bill HDDs",
      "Household FE",
      "City by month-of-sample FE",
      "City by week-of-sample FE",
      "Zip by week-of-sample FE",
      "\\textit{N}"),
    row_ends = c(
      "",
      "\\\\[-0.9em] \\midrule",
      rep("", 7)),
    sideways = F,
    title_bf = "Second-stage results",
    title_nb = "\\newline Robustness to specification: Avg. mrg. price instrumented with spot price",
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_more = NULL,
    col_numbers = T,
    col_labels = NULL,
    the_caption = paste(cap_std, cap_lag2, cap_avgmrg, cap_sig))

# tab:stage2_base_rob: Second stage, avg price, lag 2, vary HDDs and FEs -------
  tmp_mat <- lapply(X = c(
    "LogPriceBase__Lag2__CityYm_HhId__All__All__1.rds",
    "LogPriceBase__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
    "LogPriceBase__Lag2__CityYw_HhId__All__All__BillHddDays.rds",
    "LogPriceBase__Lag2__ZipYw_HhId__All__All__BillHddDays.rds",
    NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # Determine whether we have HDDs
      nc <- grepl("hdd", a_file, ignore.case = T) + 1
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[2]]$coef[nc,1] %>% cf()
      the_se <- tmp[[2]]$coef[nc,2] %>% se()
      the_cf %<>% add_stars(tmp[[2]]$coef[nc,4])
      the_f <- tmp[[1]]$f %>% n()
      the_n <- tmp[[1]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf, the_se, the_f, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Add more rows
  tmp_mat <- rbind(
    tmp_mat[1:3,],
    # HDDs
    c("F", rep("T", 3)),
    # Household FE
    c("T", "T", "T", "T"),
    # City month-of-sample FE
    c("T", "T", "F", "F"),
    # City week-of-sample FE
    c("F", "F", "T", "F"),
    # Zip week-of-sample FE
    c("F", "F", "F", "T"),
    # N
    tmp_mat[4,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage2_base_rob.tex",
    label = "tab:stage2_base_rob",
    the_mat = tmp_mat,
    row_names = c(
      "Log(Baseline price)",
      "\\quad\\textit{instrumented}",
      "First-stage F stat.",
      "Bill HDDs",
      "Household FE",
      "City by month-of-sample FE",
      "City by week-of-sample FE",
      "Zip by week-of-sample FE",
      "\\textit{N}"),
    row_ends = c(
      "",
      "\\\\[-0.9em] \\midrule",
      rep("", 7)),
    sideways = F,
    title_bf = "Second-stage results",
    title_nb = "\\newline Robustness to specification: Baseline price instrumented with spot price",
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_more = NULL,
    col_numbers = T,
    col_labels = NULL,
    the_caption = paste(cap_std, cap_lag2, cap_base, cap_sig))

# tab:stage2_sim_rob: Second stage, avg price, lag 2, vary HDDs and FEs --------
  tmp_mat <- lapply(X = c(
    "LogMrgSimIv1014__Lag2__CityYm_HhId__All__All__1.rds",
    "LogMrgSimIv1014__Lag2__CityYm_HhId__All__All__BillHddDays.rds",
    "LogMrgSimIv1014__Lag2__CityYw_HhId__All__All__BillHddDays.rds",
    "LogMrgSimIv1014__Lag2__ZipYw_HhId__All__All__BillHddDays.rds",
    NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(paste0(dir_results, a_file))
      # Determine whether we have HDDs
      nc <- grepl("hdd", a_file, ignore.case = T) + 1
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[2]]$coef[nc,1] %>% cf()
      the_se <- tmp[[2]]$coef[nc,2] %>% se()
      the_cf %<>% add_stars(tmp[[2]]$coef[nc,4])
      the_f <- tmp[[1]]$f %>% n()
      the_n <- tmp[[1]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf, the_se, the_f, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Add more rows
  tmp_mat <- rbind(
    tmp_mat[1:3,],
    # HDDs
    c("F", rep("T", 3)),
    # Household FE
    c("T", "T", "T", "T"),
    # City month-of-sample FE
    c("T", "T", "F", "F"),
    # City week-of-sample FE
    c("F", "F", "T", "F"),
    # Zip week-of-sample FE
    c("F", "F", "F", "T"),
    # N
    tmp_mat[4,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage2_sim_rob.tex",
    label = "tab:stage2_sim_rob",
    the_mat = tmp_mat,
    row_names = c(
      "Log(Simulated mrg. price)",
      "\\quad\\textit{instrumented}",
      "First-stage F stat.",
      "Bill HDDs",
      "Household FE",
      "City by month-of-sample FE",
      "City by week-of-sample FE",
      "Zip by week-of-sample FE",
      "\\textit{N}"),
    row_ends = c(
      "",
      "\\\\[-0.9em] \\midrule",
      rep("", 7)),
    sideways = F,
    title_bf = "Second-stage results",
    title_nb = "\\newline Robustness to specification: Baseline price instrumented with spot price",
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_more = NULL,
    col_numbers = T,
    col_labels = NULL,
    the_caption = paste(cap_std, cap_lag2, cap_sim, cap_sig))

# tab:cor_price: Correlations of prices ----------------------------------------
  # Load data
  tmp_mat <- readRDS(paste0(dir_results, "cor_price.rds")) %>% round(4)
  # Drop upper triangle
  for (i in 1:5) {
    for (j in 1:5) {
      if (j > i) tmp_mat[i,j] <- ""
    }
  }
  # Create the table
  simple_table(
    path = dir_output,
    filename = "cor_price.tex",
    label = "tab:cor_price",
    the_mat = tmp_mat,
    row_names = c(
      "Marginal",
      "Average",
      "Avg. Mrg.",
      "Baseline",
      "Sim. mrg.",
      NULL),
    row_ends = c(rep("\\\\[-0.5ex]", 4), ""),
    sideways = F,
    title_bf = "Price correlation",
    title_nb = "Bivariate correlations between types of prices",
    above_label = NULL,
    top_more = "  & \\multicolumn{5}{c}{\\textbf{Type of Price}} \\\\ \\cmidrule{2-6}",
    col_numbers = F,
    col_labels = c(
      "",
      "Marginal",
      "Average",
      "Avg. Mrg.",
      "Baseline",
      "Sim. mrg.",
      NULL),
    the_caption = paste(cap_avg, cap_avgmrg, cap_base, cap_sim))

# tab:stage2_neighbors: Incrementally adding neighboring zip codes -------------
  tmp_mat <- lapply(X = c(
    paste0(dir_results,  "LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds"),
    paste0(dir_results2, "Group1/LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds"),
    paste0(dir_results2, "Group2/LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds"),
    paste0(dir_results2, "Group3/LogPriceMrg__Lag2__CityYm_HhId__All__All__BillHddDays.rds"),
    NULL),
    FUN = function(a_file) {
      # Load the results
      tmp <- readRDS(a_file)
      # Determine whether we have HDDs
      nc <- grepl("hdd", a_file, ignore.case = T) + 1
      # Grab and format the coeficient and standard error
      the_cf <- tmp[[2]]$coef[nc,1] %>% cf()
      the_se <- tmp[[2]]$coef[nc,2] %>% se()
      the_cf %<>% add_stars(tmp[[2]]$coef[nc,4])
      the_f <- tmp[[1]]$f %>% n()
      the_n <- tmp[[1]]$N %>% n()
      # Return a column matrix of the results
      matrix(c(the_cf, the_se, the_f, the_n), ncol = 1)
    }) %>% do.call("cbind", .)
  # Add more rows
  tmp_mat <- rbind(
    tmp_mat[1:3,],
    # HDDs
    rep("T", 4),
    # Household FE
    rep("T", 4),
    # City month-of-sample FE
    rep("T", 4),
    # Neighbor group
    c(0:3),
    # N
    tmp_mat[4,])
  # Create the table
  simple_table(
    path = dir_output,
    filename = "stage2_neighbors.tex",
    label = "tab:stage2_neighbors",
    the_mat = tmp_mat,
    row_names = c(
      "Log(Marginal price)",
      "\\quad\\textit{instrumented}",
      "First-stage F stat.",
      "Bill HDDs",
      "Household FE",
      "City by month-of-sample FE",
      "Levels of neighboring zip codes",
      "\\textit{N}"),
    row_ends = c(
      "",
      "\\\\[-0.9em] \\midrule",
      rep("", 6)),
    sideways = F,
    title_bf = "Second-stage results",
    title_nb = "Extending the set of zip codes to neighboring zip codes",
    above_label = "\\textbf{Dependent variable:} Log(Consumption, daily avg.)",
    top_more = " & \\multicolumn{4}{c}{\\textbf{Marginal Price}} \\\\ \\cmidrule{2-5}",
    col_numbers = T,
    col_labels = c("", "Common Zips", "Neighbors 1", "Neighbors 2", "Neighbors 3"),
    the_caption = paste(cap_nbr, cap_std, cap_lag2, cap_sig))

# NOTE: I manually update this table: it was too wide to print on a page.
# # tab:balance: Balance on observables ------------------------------------------
#   # Load data
#   tmp_mat <- readRDS(paste0(dir_results, "balanceMatrix.rds"))
#   # Create the table
#   simple_table(
#     path = dir_output,
#     filename = "balance.tex",
#     label = "tab:balance",
#     the_mat = tmp_mat,
#     row_names = c(
#       "Therms consumed", "",
#       "Days in bill", "",
#       "Allowance", "",
#       "Total bill", "",
#       "HDDs", "",
#       "\\textit{N}",
#       NULL),
#     row_ends = c(rep(c("[-0.4ex]", "[1ex]"), 5), ""),
#     sideways = F,
#     landscape = T,
#     title_bf = "Balance on observables",
#     title_nb = "Comparing utilities' customers across the discontinuity",
#     above_label = NULL,
#     top_more = c(
#       "  & \\multicolumn{6}{c}{\\textbf{Summer}} & \\multicolumn{6}{c}{\\textbf{Winter}} \\\\ \\cmidrule(lr{1pt}){2-7} \\cmidrule(lr{1pt}){8-13} \n",
#       "  & \\multicolumn{3}{c}{\\textbf{Non-CARE}} & \\multicolumn{3}{c}{\\textbf{CARE}} & \\multicolumn{3}{c}{\\textbf{Non-CARE}} & \\multicolumn{3}{c}{\\textbf{CARE}} \\\\ \\cmidrule(lr{1pt}){2-4} \\cmidrule(lr{1pt}){5-7} \\cmidrule(lr{1pt}){8-10} \\cmidrule(lr{1pt}){11-13} \n",
#       NULL),
#     col_numbers = F,
#     col_labels = c(
#       "\\textbf{Variable}",
#       rep(c("PG\\&E", "SoCalGas", "\\textbf{Diff.}"), 3),
#       NULL),
#     the_caption = "Unbracketed values provide the variables' means; bracketed values denote the variables' standard deviations. The standard deviations below the difference column (\\textit{Diff.}) are pooled across utilities.")
